function y=Hessian_star(beta, Data,n)

N=max(Data(:,1));
task=max(Data(:,2));
option=max(Data(:,3));
L=size(beta,2);
x=zeros(L,L);                         %FOCs to be summed up in the loop
z=zeros(L,L);
y=zeros(L,L);
omega=.1;
b=1/N;

id=n;
for taskid=1:task
  if sum(Data((id-1)*task*option+(taskid-1)*option+1:(id-1)*task*option+taskid*option,4))~=12
    options=Data((id-1)*task*option+(taskid-1)*option+1:(id-1)*task*option+taskid*option,3);
    rank=Data((id-1)*task*option+(taskid-1)*option+1:(id-1)*task*option+taskid*option,4);
    w=Data((id-1)*task*option+(taskid-1)*option+1:(id-1)*task*option+taskid*option-1,5:5+L-option);
    w=[eye(option-1) w]; 

    %define the numerators of the likelihood
    numer(1:option-1,1)=exp(w*beta');
    %numer(4)=1;

    %define the denominator 
    denom=sum(numer)+1;

    %set up the likelihood vector
    probi=numer/denom;

    %set up the dprob/dbeta vector (option-1 by L matrix)

    %set up a vector to be used in the SOC
    for l=1:L
        expprobi(l)=w(:,l)'*numer/denom;
    end
    for op=1:option-1
       dpdb(op,:)=probi(op)*(w(op,:)-expprobi); 
    end

    %set up the SOCs
    x=x+w'*dpdb;
  end
end
x=(1-omega)*x;


for id=1:N
   for taskid=1:task
      if sum(Data((id-1)*task*option+(taskid-1)*option+1:(id-1)*task*option+taskid*option,4))~=12
        options=Data((id-1)*task*option+(taskid-1)*option+1:(id-1)*task*option+taskid*option,3);
        rank=Data((id-1)*task*option+(taskid-1)*option+1:(id-1)*task*option+taskid*option,4);
        w=Data((id-1)*task*option+(taskid-1)*option+1:(id-1)*task*option+taskid*option-1,5:5+L-option);
        w=[eye(option-1) w]; 
        
        %define the numerators of the likelihood
        numer(1:option-1,1)=exp(w*beta');
        %numer(4)=1;

        %define the denominator 
        denom=sum(numer)+1;
        
        %set up the likelihood vector
        probi=numer/denom;
        
        %set up the dprob/dbeta vector (option-1 by L matrix)
        
        %set up a vector to be used in the SOC
        for l=1:L
            expprobi(l)=w(:,l)'*numer/denom;
        end
        for op=1:option-1
           dpdb(op,:)=probi(op)*(w(op,:)-expprobi); 
        end
        
        %set up the SOCs
        z=z+w'*dpdb;
      end
   end
end

z=omega*b*z;

y=inv(x+z);